/// Atom.cpp: implementation of the atom class.
//
//////////////////////////////////////////////////////////////////////

#include "Atom.h"
#include "nr_c/nr.h"

#define PI   3.1415926535897932385

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

// what do the variables stand for?
// whats pack and position?
Atom::Atom(double InR, double InPhi, double InTheta,
		   double *pR, double *pPhi, double *pTheta, 
		   Pack   *DR, Pack   *DPhi, Pack   *DTheta,		 
		   Position *A1, Position *A2, Position *A3, Position *A4)
{
	BondLength		= InR;
	CosPhi			= cos(InPhi);
	SinPhi			= sin(InPhi);
	CosTheta		= cos(InTheta);
	SinTheta		= sin(InTheta);

	R               = pR;
	Phi				= pPhi;	
	Theta			= pTheta;

	Atom1			= A1;
	Atom2			= A2;
	Atom3			= A3;
	Atom4			= A4;

	RDerivative     = DR;
	PhiDerivative   = DPhi;
	ThetaDerivative = DTheta;

	DirectionOfPhi  = 1.0;

}


void Atom::Build(void)
{
//	cout << "Hello in atom::build now\n";
	MakeParameters();
	MakeBasis();
	MakeAtom4();
	MakeDerivativeVectors();
}

void Atom::MakeParameters(void)
{
	if(R!=NULL)
		BondLength = *R;

	if(Phi!=NULL)
		{
			
		*Phi   = fmod(*Phi, 2.0*PI);

		CosPhi = cos(*Phi);
		SinPhi = fabs(sin(*Phi), DirectionOfPhi);
		}

	if(Theta!=NULL)
		{
		*Theta   = fmod(*Theta, 2.0*PI);

		CosTheta = cos(*Theta);
		SinTheta = sin(*Theta);
		}
}


void Atom::MakeBasis(void)
{
	ZAxis = *Atom2 - *Atom3;
	ZAxis.Normalise();

	YAxis = ZAxis^(*Atom1 - *Atom2);
	YAxis.Normalise();

	XAxis = YAxis^ZAxis;

	// to understand Johns code
/*
	cout << "First atom calculated\n" << "XAxis : " << XAxis << "\n"
		<< "YAxis : " << YAxis << "\n"
		<< "ZAxis : " << ZAxis << "\n"; */
}


void Atom::MakeAtom4(void)
{
	Position Difference;

	Difference =  SinPhi*CosTheta*XAxis 
				- SinPhi*SinTheta*YAxis 
				+ CosPhi*ZAxis;

	*Atom4 = *Atom3 + BondLength*Difference;

	// to understand Johns code
/*
	cout << "In Atom::MakeAtom4\n" << "*Atom3 : " << *Atom3 << "\n"
		<< "BondLength : " << BondLength << "\n"
		<< "BondLength*Difference : " << BondLength*Difference << "\n"
		<< "*Atom4 : " << *Atom4 << "\n";
	exit(1);*/

}


void Atom::MakeDerivativeVectors(void)
{
	if(RDerivative!=NULL)
		{
		RDerivative->Vector1.Zero();

		RDerivative->Vector2 = *Atom4 - *Atom3;
		RDerivative->Vector2.Normalise();
		}

	if(PhiDerivative!=NULL)
		{
		PhiDerivative->Vector1 = (*Atom2 - *Atom3)^(*Atom4 - *Atom3);
		PhiDerivative->Vector1.Normalise(DirectionOfPhi);

		PhiDerivative->Vector2 = *Atom3^PhiDerivative->Vector1;
		}
	
	if(ThetaDerivative!=NULL)
		{
		ThetaDerivative->Vector1 = *Atom3 - *Atom2;
		ThetaDerivative->Vector1.Normalise();

		ThetaDerivative->Vector2 = *Atom3^ThetaDerivative->Vector1;
		}
}

